Experiment_Preprocessing

to do before publish

#change all path to data showing my environment (NAS...) mettre lien à partir de Workspace ou il dowloaderont les fichier source (fastq et objet R)
#ajouter un test pour etre sur que y'ai subset 1 et 2 et pas que un des deux 

#do we want to show the code in rmd ? 

# GO mise en ligne des fichier brut

#test final sans aucun objet, avec le premier obj et avec le deuxieme

#change output path

[1] “You are starting the analysis from count matrix obtain with CellRanger and Cite-seq-count”

Loading the first experiment

Looking to HTO distribution accross sample :

#HTO distribution
par(mfrow=c(1,1))
par(las=2)
par(mar=c(5,15,3,3))
barplot(rowSums(hto1), main = "sequenced HTO distribution", horiz=TRUE)

rowSums(hto1)

##   Spleen-MP    Spleen-M Spleen-ctrl    Spleen-P   Thymus-MP    Thymus-M 
##      688432      206663      265412      662165      196262      117402 
## Thymus-ctrl    Thymus-P 
##      156856      562594

Demultiplexing HTO

Demultiplexing results

Cells classification

print (table(hashtag1@meta.data$MULTI_ID))

## 
##     Doublet    Negative Spleen-ctrl    Spleen-M   Spleen-MP    Spleen-P 
##         307          69         438         255         710         609 
## Thymus-ctrl    Thymus-M   Thymus-MP    Thymus-P 
##         325         245         294         935

Violinplot (features)

VlnPlot(hashtag1,features = c("nFeature_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

Violinplot (HTO counts)

VlnPlot(hashtag1,features = c("nCount_HTO"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

tSNEs based on HTO

# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = hashtag1, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
hashtag1 <- RunTSNE(hashtag1, distance.matrix = hto.dist.mtx, perplexity = 100)

HTO margin

Tsne<-data.frame(
  tSNE_1 = hashtag1@reductions$tsne@cell.embeddings[,1],
  tSNE_2= hashtag1@reductions$tsne@cell.embeddings[,2],
  gene= hashtag1@meta.data$HTO_margin
)

HTO= hashtag1@meta.data$MULTI_ID
Max=max(hashtag1@meta.data$HTO_margin)
Min=min(hashtag1@meta.data$HTO_margin)
ggplot(Tsne,aes(x=tSNE_1,y=tSNE_2))+geom_point(aes(color=gene,shape=HTO))+
           scale_colour_gradient2(low = "blue",mid="orange",high="red",name="HTO margin",midpoint=(Max+Min)/2)+scale_shape_manual(values = c(15,16,17,18,19,20,21,22,23,24,25))

Ridge plots

Visualize enrichment for selected HTOs with ridge plots

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[1:2],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[3:4],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[5:6],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag1, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[7:8],ncol = 2, group.by = "MULTI_ID")

Sample Information

The analysis will be run on the sample 1 (181031).

During the sample loading, we filter cells that do not pass the following filters.
Here are the description of those parameters in the Seurat CreateSeuratObject function:

• min.genes: Include cells where at least 200 genes are detected
• min.cells: Include genes with detected expression in at least 3 cells

After those filters, the remaining cell number is `r length(colnames(hashtag1@assays$RNA@data)).

#add Exp1 cell identity
HTO_cr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-ctrl" ))
HTO_ct1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-ctrl" ))
HTO_mr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-M" ))
HTO_mt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-M" ))
HTO_pr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-P" ))
HTO_pt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-P" ))
HTO_pmr1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Spleen-MP" ))
HTO_pmt1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Thymus-MP" ))
HTO_d1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Doublet" ))
HTO_n1 <- row.names(subset(hashtag1@meta.data, MULTI_ID == "Negative" ))

HTO_thymus1 = c(HTO_ct1,HTO_mt1,HTO_pt1,HTO_pmt1)
HTO_spleen1 = c(HTO_cr1,HTO_mr1,HTO_pr1,HTO_pmr1)
HTO_identified1 = c(HTO_thymus1, HTO_spleen1)


# Create a Seurat object without doublet and unassigned cells (remove "negative", "doublet" & "nothing"))
clean.subset1 <- subset(x = hashtag1, cells = HTO_identified1)
VlnPlot(clean.subset1,features = c("nFeature_RNA", "nCount_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

After removing doublets and negative cells, the remaining cell number is `r length(colnames(clean.subset1@assays$RNA@data)).

Mitochondrial percentage versus nFeatures

df<-data.frame(multi.id=Seurat1@misc$old_meta_data$MULTI_ID,percent.mito=Seurat1@misc$old_meta_data$percent.mito,nFeature_RNA=Seurat1@misc$old_meta_data$nFeature_RNA)
ggplotly(ggplot(df,aes(x=nFeature_RNA,y=percent.mito,color=multi.id))+geom_point())

UMAP

ggplotly(DimPlot(Seurat1, reduction = "umap", group.by = "MULTI_ID", do.label = TRUE, pt.size = 1))

Loading the second experiment

par(mfrow=c(1,1))
par(las=2)
par(mar=c(5,15,3,3))
barplot(rowSums(hto2), main = "sequenced HTO distribution", horiz=TRUE)

rowSums(hto2)

##   Spleen-MP    Spleen-M Spleen-ctrl    Spleen-P   Thymus-MP    Thymus-M 
##     3569690     1668124      929963     3251282      484276      749822 
## Thymus-ctrl    Thymus-P 
##      310711     1203239

Demultiplexing results

Cells classification

print (table(hashtag2@meta.data$MULTI_ID))

## 
##     Doublet    Negative Spleen-ctrl    Spleen-M   Spleen-MP    Spleen-P 
##        2254         427         824        1294        2589        1689 
## Thymus-ctrl    Thymus-M   Thymus-MP    Thymus-P 
##         461        1266         522        1664

Violinplot (features)

VlnPlot(hashtag2,features = c("nFeature_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

Violin plots (HTO counts)

VlnPlot(hashtag2,features = c("nCount_HTO"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

tSNEs based on HTO

# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = hashtag2, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
hashtag2 <- RunTSNE(hashtag2, distance.matrix = hto.dist.mtx, perplexity = 100)

DimPlot(hashtag2, group.by = "MULTI_ID")

HTO margin

Tsne<-data.frame(
  tSNE_1 = hashtag2@reductions$tsne@cell.embeddings[,1],
  tSNE_2= hashtag2@reductions$tsne@cell.embeddings[,2],
  gene= hashtag2@meta.data$HTO_margin
)

HTO= hashtag2@meta.data$MULTI_ID
Max=max(hashtag2@meta.data$HTO_margin)
Min=min(hashtag2@meta.data$HTO_margin)

ggplot(Tsne,aes(x=tSNE_1,y=tSNE_2))+geom_point(aes(color=gene,shape=HTO))+
           scale_colour_gradient2(low = "blue",mid="orange",high="red",name="HTO margin",midpoint=(Max+Min)/2)+scale_shape_manual(values = c(15,16,17,18,19,20,21,22,23,24,25))

Ridge plots

Visualize enrichment for selected HTOs with ridge plots

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[1:2],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[3:4],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[5:6],ncol = 2, group.by = "MULTI_ID")

RidgePlot(hashtag2, assay = "HTO", features = rownames(GetAssayData(hashtag1,assay = "HTO"))[7:8],ncol = 2, group.by = "MULTI_ID")

Sample Information

The analysis will be run on the sample 2 (190211).

During the sample loading, we filter cells that do not pass the following filters.

Used parameters in the Seurat CreateSeuratObject function: * min.genes: 3 . Include cells where at least 3 genes are detected * min.cells: 200 . Include genes with detected expression in at least 200 cells

After those filters, the remaining cell number is `r length(colnames(hashtag2@assays$RNA@data)).

#add Exp2 cell identity
HTO_cr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID == "Spleen-ctrl" ))
HTO_ct2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-ctrl" ))
HTO_mr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-M" ))
HTO_mt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-M" ))
HTO_pr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-P" ))
HTO_pt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-P" ))
HTO_pmr2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Spleen-MP" ))
HTO_pmt2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Thymus-MP" ))
HTO_d2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Doublet" ))
HTO_n2 <- row.names(subset(hashtag2@meta.data, MULTI_ID== "Negative" ))

HTO_thymus2 = c(HTO_ct2,HTO_mt2,HTO_pt2,HTO_pmt2)
HTO_spleen2 = c(HTO_cr2,HTO_mr2,HTO_pr2,HTO_pmr2)
HTO_identified2 = c(HTO_thymus2, HTO_spleen2)

# Create a Seurat object without doublet and unassigned cells (remove "negative", "doublet" & "nothing"))
clean.subset2 <- subset(x = hashtag2, cells = HTO_identified2)
VlnPlot(clean.subset2,features = c("nFeature_RNA", "nCount_RNA"),pt.size = 0.1, log = TRUE,  group.by = "MULTI_ID")

After removing doublets and negative cells, the remaining cell number is `r length(colnames(clean.subset2@assays$RNA@data)).

Adding ADT

# Load in the UMI matrix
umi <- GetAssayData(object = clean.subset2, slot = "counts")

# Load in the ADT count matrix
raw.adt <- Read10X(PATH_ADT_DATA2, gene.column = 1)
adt <- raw.adt[c(1:6),]

rownames(adt) <- c("CD4","CD5","CD8","CD25","CD44","CD69")

#create an empty matrix containing NAs
Cell.list <- colnames(GetAssayData(object = clean.subset2[["RNA"]], slot = "data" ) )
ADT.list <- c(unique(rownames(adt)))
mat.adt <- matrix(nrow = length(ADT.list), ncol = length(Cell.list))
rownames(mat.adt) = ADT.list
colnames(mat.adt) = Cell.list

# Get cell barcodes detected by both RNA and ADT
joint_bcs <- intersect(colnames(umi),colnames(adt))
adt <- as.matrix(adt[,joint_bcs])

# Fill the empty matrix with values when existing
mat.adt[,joint_bcs]<-adt[,joint_bcs]

# Add ADT data as a new assay independent from RNA
clean.subset2[["ADT"]] <- CreateAssayObject(counts = mat.adt[,colnames(clean.subset2)])

# Normalize ADT data, here we use centered log-ratio (CLR) transformation
clean.subset2 <- NormalizeData(clean.subset2, assay = "ADT", normalization.method = "CLR")

#Scale
clean.subset2 <- ScaleData(clean.subset2, assay = "ADT")

ADT list :

print (rownames(adt))

[1] “CD4” “CD5” “CD8” “CD25” “CD44” “CD69”

Mitochondrial percentage versus nFeatures

df<-data.frame(multi.id=Seurat2@misc$old_meta_data$MULTI_ID,percent.mito=Seurat2@misc$old_meta_data$percent.mito,nFeature_RNA=Seurat2@misc$old_meta_data$nFeature_RNA)
ggplotly(ggplot(df,aes(x=nFeature_RNA,y=percent.mito,color=multi.id))+geom_point())

UMAP:

ggplotly(DimPlot(Seurat2, reduction = "umap", group.by = "MULTI_ID", do.label = TRUE, pt.size = 1))

Merging our two experiments

Load separate R object

You can load objects done with the code above. Or our object ?? (link )

Integrating the 2 seurat objects with seurat integration (cca)

We identified r length(gene1) expressed in sample1 and r length(gene2) expressed in sample2. r length(common_genes) are in common in this two set and will the integrated in the merged and corrected object.

UMAP:

Analysis part

Sample Information

The analysis will be run on the sample MYC_PTEN_01 (181031 & 190211).

During the sample loading, we filter cells that do not pass the following filters. We also filter cells that are detected as human/mouse multiplet using their barcodes.
Here are the description of those parameters in the Seurat CreateSeuratObject function:

• min.genes: Include cells where at least this many genes are detected
• min.cells: Include genes with detected expression in at least this many cells

After those filters, and merging MYC_PTEN_01 and MYC_PTEN02 the remaining cell number is r length(colnames(exp1.2.integrated)).

UMAP and clustering parameter

Merge checking

HTO

Orig.idents

T-cell selection

According to T-cell markers we will exclude Cd3d low clusters: 13 (Bcells), 11, 14, 17 (monocytes/macrophages). According to T-cell markers we will exclude Cd3d/Cd3e low clusters: 11 (Bcells), 13, 17, 18, 19 (monocytes/macrophages), 16, 14 (ILC/NK).

Known RNA B and T markers

re clustering

T-cell umaps

HTO

clustering

END OF PREPROCESSING

We obtain the final object with clustering to start the analysis

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.18.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.23         RColorBrewer_1.1-2 magrittr_1.5      
## [4] dplyr_0.8.1        gridExtra_2.3      kableExtra_1.1.0  
## [7] plotly_4.9.0       ggplot2_3.1.1      Seurat_3.0.1      
## 
## loaded via a namespace (and not attached):
##   [1] nlme_3.1-140        tsne_0.1-3          bitops_1.0-6       
##   [4] webshot_0.5.1       httr_1.4.0          sctransform_0.2.0  
##   [7] tools_3.5.3         R6_2.4.0            irlba_2.3.3        
##  [10] KernSmooth_2.23-15  lazyeval_0.2.2      colorspace_1.4-1   
##  [13] npsurv_0.4-0        withr_2.1.2         tidyselect_0.2.5   
##  [16] compiler_3.5.3      rvest_0.3.4         xml2_1.2.0         
##  [19] labeling_0.3        caTools_1.17.1.2    scales_1.0.0       
##  [22] lmtest_0.9-37       readr_1.3.1         ggridges_0.5.1     
##  [25] pbapply_1.4-0       stringr_1.4.0       digest_0.6.19      
##  [28] rmarkdown_1.12      R.utils_2.8.0       base64enc_0.1-3    
##  [31] pkgconfig_2.0.2     htmltools_0.3.6     bibtex_0.4.2       
##  [34] highr_0.8           htmlwidgets_1.3     rlang_0.3.4        
##  [37] rstudioapi_0.10     shiny_1.3.2         zoo_1.8-5          
##  [40] jsonlite_1.6        crosstalk_1.0.0     ica_1.0-2          
##  [43] gtools_3.8.1        R.oo_1.22.0         Matrix_1.2-17      
##  [46] Rcpp_1.0.1          munsell_0.5.0       ape_5.3            
##  [49] reticulate_1.12     R.methodsS3_1.7.1   stringi_1.4.3      
##  [52] yaml_2.2.0          gbRd_0.4-11         MASS_7.3-51.4      
##  [55] gplots_3.0.1.1      Rtsne_0.15          plyr_1.8.4         
##  [58] grid_3.5.3          promises_1.0.1      parallel_3.5.3     
##  [61] gdata_2.18.0        listenv_0.7.0       ggrepel_0.8.1      
##  [64] crayon_1.3.4        lattice_0.20-38     cowplot_0.9.4      
##  [67] splines_3.5.3       hms_0.4.2           SDMTools_1.1-221.1 
##  [70] pillar_1.4.0        igraph_1.2.4.1      future.apply_1.2.0 
##  [73] reshape2_1.4.3      codetools_0.2-16    glue_1.3.1         
##  [76] evaluate_0.13       lsei_1.2-0          metap_1.1          
##  [79] data.table_1.12.2   httpuv_1.5.1        png_0.1-7          
##  [82] Rdpack_0.11-0       gtable_0.3.0        RANN_2.6.1         
##  [85] purrr_0.3.2         tidyr_0.8.3         future_1.13.0      
##  [88] assertthat_0.2.1    xfun_0.7            mime_0.6           
##  [91] rsvd_1.0.0          xtable_1.8-4        later_0.8.0        
##  [94] survival_2.44-1.1   viridisLite_0.3.0   tibble_2.1.1       
##  [97] cluster_2.0.9       globals_0.12.4      fitdistrplus_1.0-14
## [100] ROCR_1.0-7